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Abstract. - We demonstrate that the quantum dynamics of a many-body Fermi-Bose system can 
be simulated using a Gaussian phase-space representation method. In particular, we consider the 
apphcation of the mixed fermion-boson model to ultracold quantum gases and simulate the dynam- 
ics of dissociation of a Bose-Einstein condensate of bosonic dimers into pairs of fermionic atoms. 
We quantify deviations of atom-atom pair correlations from Wick's factorization scheme, and show 
that atom-molecule and molecule-molecule correlations grow with time, in clear departures from 
pairing mean-fleld theories. As a first-principles approach, the method provides benchmarking of 
approximate approaches and can be used to validate dynamical probes for characterizing strongly 
correlated phases of fermionic systems. 



C/3 



> 
o 



On 
O 



X 



. Introduction. — The physics of interacting fermions 
ife the basis of many of the most important phenomena in 
condensed matter physics, ultracold gases, and quantum 
chemistry. A fundamental issue is how the microscopic in- 
teractions at the quantum level give rise to collective and 
Emergent effects in many-body systems. For many situa- 
tions, particularly in condensed matter systems, static or 
Equilibrium correlation functions are sufficient to connect 
theory and experiment, and sophisticated techniques have 
been developed to calculate and measure them [T]. 

' Addressing similar questions in the domain of many- 
body dynamics, however, has limitations in condensed 
matter systems. Ultracold quantum gases, on the other 
Ijand, allow creation of highly controllable implementa- 
tions of analogue many-body systems for which the dy- 
namical evolution and correlations are directly accessible 
PHB]. The purity and tunability of these 'tailor-made' 
analogue systems means that ultracold quantum gases are 
ideal for testing fundamental ideas in quantum many-body 
physics and are leading candidates for dynamical 'quan- 
tum simulation' [THlOj . In order to make predictions from 
the underlying theory and to validate the potential sim- 
ulators [nHH], or to benchmark approximate approaches, 
a numerical simulation of the real-time dynamics is re- 



quired. Similar requirements of exact simulation of many 
fermions arise in determining the quantum chemistry of 
complex molecular systems [12] . 

In this work we perform first-principles dynamical sim- 
ulations of a fermion-boson model. We use a Gaussian 
stochastic method based on a generalized phase-space rep- 
resentation of the quantum density operator |13j . The 
fermion-boson model forms the underlying basis for a 
broad range of phenomena in condensed matter and ul- 
tracold atom physics. It was originally proposed in the 
context of high-temperature superconductivity [T3], but 
in ultracold gases it corresponds to the theory of reso- 
nance superfiuidity with Feshbach molecules [15] . The lat- 
ter forms the basis of a two-channel model for describing 
the physics of the BCS-BEC crossover [TB]. More recently, 
the fermion-boson model has been used for analyzing the 
decay of double occupancies (doublons) |17| in a driven 
Fermi- Hubbard system |18j . The particular situation that 
we simulate here corresponds to spontaneous dissociation 
of a Bose-Einstein condensate (BEG) of molecular dimers 
into fermionic atoms P^[^ , in which case the model pro- 
vides the fermionic equivalent of parametric downconver- 
sion in quantum optics: the production of pairs of entan- 
gled particles. 
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The Gaussian phase-space method represents an exten- 
sion to fermions of successful bosonic techniques pTHSO] . 
The essence of the method is the mapping of the density 
operator evolution onto a Fokker-Planck equation for a 
phase-space distribution, via a continuous Gaussian op- 
erator basis |13| . The evolving distribution is then sam- 
pled with stochastic differential equations (SDEs) for the 
phase-space variables [1T1[22]- The mapping to the phase- 
space distribution is exact [21"H23ll3Tj provided no bound- 
ary terms arise in deriving the the Fokker-Planck equation. 
In practice, such terms may develop after some simulation 
time, but in an easily verifiable way [131I3T], putting a 
well characterised upper limit to a useful simulation du- 
ration. Numerical signatures of systematic errors include: 
[i) the onset of spiking behaviour due to the presence of 
near-singular trajectories; [ii) sudden dramatic increase 
in the statistical uncertainties; and {in) development of 
power-law tails in the probability distribution. All these 
signatures, well-known from the early studies of real-time 
dynamics of bosonic systems [23] and from equilibrium cal- 
culations for fermion systems using imaginary-time tech- 
niques [131131] . carry over to the present simulations of 
fermion dynamics and are verified in the numerical exam- 
ples that we present below. 

As in any stochastic method, sampling error limits the 
precision of the results. However, unless the distribu- 
tion develops power-law tails, indicated by the above- 
mentioned signatures, this uncertainty can be made ar- 
bitrarily small by increasing the number of trajectories. 

From the physical point of view, the Gaussian phase- 
space method can be viewed as providing the quantum 
corrections, through additional stochastic terms, to differ- 
ent mean- field approaches. For example, with certain fac- 
torization assumptions I32| , the method reduces to a time- 
dependent Hartree-Fock formalism. Furthermore, neglect- 
ing the stochastic terms recovers the approximate pairing 
mean-field theory (PMFT) [19j[33], to which we compare 
the phase-space results. While often accurate for deter- 
mining particle number densities, the mean- field approach 
gives no direct information about higher-order correla- 
tions, and its accuracy is not known a priori. In contrast 
to this, the first-principles simulations presented here re- 
veal significant development of higher-order correlations. 

For the first application of the fermionic phase-space 
method to a multimode dynamical problem, we consider 
a uniform molecular BEG (MBEC) initially in a coherent 
state at zero temperature, with no atoms present. Assum- 
ing sufficiently low densities, we neglect s-wave scattering 
interactions to simplify the treatment. The Hamiltonian 
of this fermion-boson model [14] is given by 



a^mk 



rhta ] , (1) 



where k labels the plane- wave modes and cr = 1,2 labels 
the effective spin state for the atoms. Even though we 
will present the numerical results for a one-dimensional 
(ID) system, we formulate the problem in the general 



case as the method is straightforward to use in higher 
dimensions. The fermionic number and pair operators 
are defined as nk,cr = cj. ^Ck.o- and Wk = Ck,ic_k,2, with 
{ck,cr, cj^/ g.' }='^kk"5CTCT' , while the bosonic molecular opera- 
tor obeys [a, a^] = 1. The atom- molecule coupling (invoked 
by a magnetic Feshbach resonance sweep or optical Raman 
transitions) is characterized by k = xd/L^^"^ [331, where 
L is the size of the quantization box, and mediates an 
effective interaction between the atoms. The first term, 
hlS.\^ = h? |k|^ I (2m a) -I- ?iA, contains the kinetic energy of 
the atoms (of mass ma), while the detuning A < cor- 
responds to the total dissociation energy 2/i|A| imparted 
onto the system by the external fields. 

Because of the symmetry between spins in the Hamil- 
tonian, and the equal initial populations, we need only to 
consider the number operator for one of the spin states 
Tik = fi-k = "-k.i = ^k,2- An additional operator identity 
that follows from the Hamiltonian is 



TOfcTOk (= fik,ifi-k,2) = fik, 



(2) 



which arises because the condensate to which the atom 
pairs arc coupled is assumed to be homogeneous. One 
consequence of eq. (I2|) is that the relative number of 
atoms with equal and opposite momenta is perfectly 
squeezed |20j . i.e. with zero variance. It also means 
that the second-order atom-atom correlation function re- 
duces to (k, -k) = (mj^mk)/(fik,i)("-k,2) = l/("-k)- 
Thus the atom-atom correlation function can be deter- 
mined from the number density alone. 

One effective approximate approach for treating the dy- 
namics of dissociation is the PMFT [THIE^ , which is ob- 
tained by assuming atom-molecule decorrelation and by 
replacing the molecular operator by a coherent mean-field 
amplitude, a — > /3. 

In this paper we solve the full Hamiltonian ([1]) exactly, 
and in order to quantify deviations from the PMFT be- 
havior we evaluate several correlation functions. The de- 
partures from Wick decorrelation are analyzed via the cor- 
relation coefficient 



(3) 



which is unity within the PMFT. 

To examine molecule-atom pair correlations and the 
second-order coherence of the molecular field, we define 



5L'i(k) 



a^dnk) 



a 



?^k; 



(a^ d)aa) 



(4) 



Again, within the PMFT, these will be unity. We may 
expect that, over time, correlations will develop between 
the molecular and atomic fields; the Gaussian phase-space 
simulations give exact quantitative accounts of these ef- 
fects. 

Fermionic phase-space representation. As 

mentioned above, the essence of the Gaussian phase- 
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space method is the mapping of the density opera- 
tor evolution into a set of stochastic differential equa- 
tions [12 that can be solved on a computer. For 
a large class of second-quantized system Hamiltonians, 
such as those containing no higher than quartic terms 
in field operators, such mapping introduces no addi- 
tional approximations. For M plane-wave modes, the 
quantum state governed by Hamiltonian ([l} can be 
mapped onto a complex phase-space of dimension 3il/ + 
2: A(<) = (ni,...,riM,TOi,.--,mA/,TO;^,...,m|^,/3,/3+), 
with ^ m* and (3'^ ^ (3*. The phase-space equations 
are structurally similar to the Heisenberg equations for the 
corresponding operators, but the stochastic terms are not 
unique and can be modified by different choices of diffu- 
sion gauge [T3ll2ll[34] . This freedom allows the equations 
to be tailored to have different numerical properties. One 
specific set of Ito SDEs is |M] : 

nk = am^ + a+mu + ^/^^^k ("ikCi + "^k C2) 

rhk = -2i(5kmk + a (1 - 27ik) + (^kCi - ^iQ) 

m+ = 2*<5km+ + a+ (f - 2nk) + ^ {m+\i - nld) 

« = -i^Ek"^k + 7^Ci 
"+ = -w;7Ek"i^ + 7^C2, 

(5) 

where the derivative is with respect to a scaled time, 
r = t/^O: with to = 1/k-\//Vo. We have normalized 
the molecular field by its maximum (initial) value, a = 
P/V^, where Nq is the initial number of molecules. 
The complex Gaussian noises Q obey {(j{T)(j' {t')) = 
0, (0(t)C*'(t')) = Sjf6{T - t'). This form of eqs. © 
shows that with drift terms of order 1, the noise terms 
are ^ Nq , i.e. the noise and therefore non-mean- 
ficld corrections to correlations become more important 
for decreasing A^o- In practice we convert the equations 
to Stratonovich form and integrate with a semi-implicit 
method. Stochastic averages of the variables give the first- 
order operator moments; normally ordered higher-order 
moments are obtained by averages of the corresponding 
Wick decomposition [T3|, e.g. {ml^rhk) = {m^m\^)s + 
{n^)s- Note that the final, averaged moment will not sat- 
isfy Wick's theorem for a general quantum state. 

The stochastic sampling assumes a sufficiently bounded 
distributions, such that any boundary terms could have 
been neglected in obtaining the Fokker-Planck equation. 
Previous experience with bosons [22 and fermions [3T] 
has shown that spikes and a sudden rapid growth of the 
statistical sampling errors in observables are seen when 
the tails of the probability distributions do not decay fast 
enough. Through comparisons with independent numer- 
ical solutions we confirm that the onset of such spiking 
behaviour and the rapidly growing sampling error signifies 
the limit of the useful simulation time. Although it can 
be controlled somewhat with gauges, [T31IM1I31] , the finite 
simulation time is a well-known limitation of stochastic 
phase-space methods. In the numerical examples of figs. 
(HO and m below, the simulation results are shown for time 



windows well below the spiking time. We simulate a suf- 
ficiently large number of stochastic trajectories to reduce 
the statistical errors to below the thickness of the relevant 
curves shown in the figures. In addition, we use identities 
such as ^ and different gauges as a further demonstra- 
tion that the simulations are exact up until the emergence 
of spiking behaviour. 

Unlike quantum Monte Carlo approaches that are well- 
suited to calculation of exact ground-state properties or 
to simulation through imaginary time j35| . the Gaussian 
stochastic method does not suffer from a 'dynamical sign 
problem' |36j . Other approaches for real-time simula- 
tions include the time-dependent density functional the- 
ory [37] . although in practice this method is often re- 
stricted in accuracy by the need for exact functionals. 
Methods that use matrix-product-states based algorithms 
have been very successful for applications to one spatial 
dimension j55H42| . however, as these methods require a 
truncated basis they do not fulfill the strict benchmark- 
ing criteria that a first-principles method can provide. An 
interesting direction in recent years has been the exten- 
sion to fermionic systems of stochastic wavefunction ap- 
proaches [13] J which are similar in spirit to phase-space 
methods. 

Few-mode system. — To confirm the validity of our 
numerical implementation of the phase-space method, we 
first independently solve a small system with A'o = 10 
molecules and M = 10 atomic modes in a standard 
number-state basis. For this test system, with a bosonic 
number-basis truncation of rimax ^ 10^, the Hilbert space 
has dimension d ~ 2*^?iinax — 10^. In fig. [T] we show 
the population in the momentum modes (fik) calculated 
using the phase-space method (with f 0^ stochastic trajec- 
tories) and the number-state basis; we also illustrate the 
identity ^ by calculating and plotting {rhj^^rhkg) directly 
and comparing it with {hko)- The agreement between the 
two methods is excellent. The top two curves in fig. [T] 
illustrate the deviation of the PMFT prediction (dashed- 
dotted curve) from the exact calculation (solid curve) for 
the resonant mode fco- 

In fig. [5J we show the population dynamics of the res- 
onant mode, TT-fco, together with the explicit behaviour of 
the statistical errors due to the stochastic sampling. The 
number of stochastic trajectories in this example (which 
otherwise is the same as in fig. [T]) is small enough to ren- 
der the sampling errors visible. At the same time, we have 
chosen the time window longer than the onset of spiking 
behaviour and the sudden dramatic growth of statistical 
errors to explicitly illustrate the signatures of systematic 
errors that limit the useful simulation time of the phase- 
space method. 

To further evaluate the differences between treating 
the Hamiltonian ([T]) exactly and using the approximate 
PMFT, we plot in fig. [3] (a) the correlation coefficient 
W . Clear deviations are seen as time evolves; the devi- 
ations illustrate that in the exact treatment the following 
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Fig. 1: (Color online) Comparison between different meth- 
ods for the population dynamics of individual atomic modes 
(nfe(T)): phase-space method with 10^ stochastic trajecto- 
ries (solid black curves), number-state basis (dashed yellow 
curves). The top curve with the largest oscillation ampli- 
tude is for the resonant mode ka — 6dk (corresponding to 
Sko = to^ko = 0); the other curves are the sidebands stepped 
by dk = 160~^''*/do, where the lengthscale is do = -^/ftfo/Sma. 
For the fco-mode, we also plot the result from the PMFT for 
comparison (top dash-dotted curve). To illustrate the iden- 
tity ([2]) we plot the quantity {m^^nikg) as large black squares 
for the phase-space method and as small yellow squares for 
the number-state calculation. The left inset shows the num- 
ber of molecules A''^ = {a^ a) (top curve) and the total num- 
ber of atoms Na ~ "^ki^i^) '^'i^ °f ^Pi'^ states. For all 
curves from the phase-space method, the statistical errors of 
±1(7 standard deviation are smaller than the thickness of the 
curves, within the time window shown. The inset on the right 
shows these errors explicitly as dark and light grey shadings 
corresponding to quantities {nfco(r)) and {mj.^rhko) , respec- 
tively. 

inequality holds, (mj(.TOk) > |(?TT-k)|^ + ('^k)'^, whereas the 
PMFT prescribes an equality sign. Next, we consider the 
molecule-atom and molecule-molecule second-order corre- 
lations, gina and gmm [see figs. [3] (b) and (c)]. Within 
the PMFT, both correlations are identically equal to 1. 
However, our exact results show that the molecule-atom 
correlation initially grows with time while the total atomic 
population grows. Then it changes to anti-correlation as 
the resonant-mode atoms start to re-associate. Meanwhile 
the effect on the molecular field from the atom interactions 
is revealed as it gradually loses its second-order coherence, 
albeit not by a significant amount for this few-mode sys- 
tem. 

Multi-mode system. — We now use the phase-space 
method for simulating large ID systems, with AI — 10^ 
atomic modes and Nq = 10^ — 10'' {'^^K2) molecules at den- 
sities niD — 1-3 X 10^ — 1.3 X lO'' m~^. In these cases, the 
number-state calculation is impossible as the dimension of 
the Hilbert space is enormous (o? = 2*^nmax ^ 10'^'^"). In 
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Fig. 2: (Color online) (a) Population of the resonant mode, 
(wfef, ), as a function of time r for the same parameters as in 
fig. [T] except with 10^ stochastic trajectories and a longer 
simulation time. The solid line is from the phase-space method, 
with the grey shading representing statistical errors of zblcr 
standard deviation. The dashed yellow line is from the number- 
state calculation, (b) Evolution of the statistical error on rikg , 
showing the emergence of spiking behaviour past r > 3 and 
the subsequent rapid growths of the sampling errors. The size 
of the statistical errors - prior to the emergence of the spiking 
behaviour - scales as l/V Ns as expected from the central limit 
theorem, where Ns is the number of stochastic trajectories. 
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Fig. 3: (Color online) (a), (b), (c) Correlation coefficients W, 
3ma(fco), and gmln as a function of time for the test system 
of fig. [TJ with A'^o = 10 and M = 10. Solid black curves are 
from the phase-space method; dashed yellow curves are from 
the number-state calculation. The lower curve in (a) is for 
the full summation in eq. Q, whereas the upper curve is the 
respective correlation coefficient for the resonant mode ko- (d), 
(e), (f) The same as in the left column, but from the phase- 
space method for A'o = 10^ and M = 10^. The dashed grey 
curve in (f) is from an ensemble of PMFT calculations with 
a set of Poissonian- weighted Nqj and A'^o ~ 10^ (see text). 
Statistical errors of ±lcr standard deviation for all curves from 
the phase-space method are contained within the thickness of 
the curves and are obtained from 10^ stochastic trajectories. 
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Fig. 4: (Color online) The fraction of remaining molecules, 
N{t)/No, as a function of time, for: A^o = 10**, 5 — —2.5 (top 
curve); A^o — 10^, 5 = — 25 (intermediate curve); A^'o = 10^, 
S — —2.5 (bottom curve). In all cases the coupling k is chosen 
to result in the same timescale to = 1 / k V A^o = 2 x 10"* s. The 
solid curves are from the phase-space method (with 10^ trajec- 
tories), whereas the dash-dotted curves are from the approxi- 
mate PMFT method; the difference between the two types of 
curves (seen in the inset, for the lowest curve) is almost invisi- 
ble on the scale of the graph. The curve for 5 = —25 illustrates 
a route away from the regime of strong molecular depletion by 
increasing the dissociation energy 2ft |A| by an order of magni- 
tude, while keeping the same A^o as for the lowest curve. For 
all curves from the phase-space method, the statistical errors 
of ±1(7 standard deviation are smaller than the thickness of 
the curves, within the time window shown. 



fig.Hlwc show the evolution of the number of molecules for 
three different cases. For the top curve, the initial number 
(A^o = 10'') is much larger than the number of available 
atomic modes, each of which hosts at most 1 atom due to 
the Pauli blocking. Accordingly, we see negligible deple- 
tion of the MBEC, which makes the relative size of the 
bosonic fluctuations very small. Hence, we do not observe 
significant deviations from the PMFT, including in the 
molecular second-order coherence, eq. (j?]), which differed 
from 1 by less than 10^'' in this case. 

The situation changes for the bottom curve, for which 
A^o — 10^ is comparable with the number of atomic modes 
within the relevant width of the momentum distribution 
near fcg; we estimate this number [33] to be ~ O.IM = 10^. 
In this case, we see strong molecular depletion and an in- 
creased role of bosonic quantum fiuctuations so that the 
PMFT starts to show disagreement with the exact result. 
Admittedly the disagreement is still very small, implying 
that the predictions of the PMFT for total particle num- 
bers can be rather accurate. The same is not true, how- 
ever, for higher-order correlations, shown in the right col- 
umn of fig. 13] for the same parameters as the lowest curve 
of fig. S] Here, the large depletion of the MBEC and the 
increased role of quantum fiuctuations are manifest - be- 
yond the predictability of the PMFT - in strong higher- 
order correlations. The correlation coefficient W clearly 
deviates from one, though to a lesser extend than in the 
few-mode system. The deviations of the molecule-atom 
and molecule- molecule correlations from gma{kn) = 1 and 

(2) 

gmrn = 1, On the Other hand, are more dramatic. 



The development of decoherence in the molecular field 
can largely be accounted for by the dephasing of atomic- 
molecular oscillations due to total number uncertainty. 
The frequency of the oscillations depends on the initial 
number of molecules; with a range of frequencies, the os- 
cillations get out of phase and thus, for example, prevent 
complete disassociation of the molecular field from being 
seen in the average. As illustrated by the dashed grey 
curve in fig. [3] (f), this effect can be reproduced by an 
ensemble of mean-field trajectories with different initial 
numbers. Here we have used a Poissonian weighting to 
give the same number-distribution as the initial coherent 
molecular field in the exact calculation, with A^o = 100 
[44] . The large values of gmln, which are possible with 
a state containing a superposition (or mixture) of a few 
low-occupation number states [45] , occur at maximum de- 
pletion when the number-uncertainty is relatively large. 

Summary. — We have demonstrated a successful ap- 
plication of a fermionic phase-space representation to 
first-principles quantum dynamics of a fermion-boson 
model. We simulated the coherent molecular dissocia- 
tion to fermionic atoms and found significant higher-order 
correlations that cannot be accounted for by the approx- 
imate pairing mean-field theory. The knowledge of such 
correlations and the development of experimental probes 
to measure them provide the most accurate characteriza- 
tion of quantum many-body phases in strongly correlated 
systems. 

The method is exact up until clear 'spiking' signatures 
emerge in the stochastic trajectories. For the present 
model, the useful simulation duration encompasses the 
time required for partial recombination of the atoms into 
molecules and for significant higher-order correlations to 
emerge. The accuracy during this useful simulation time 
was independently confirmed by comparison to number 
state calculations (for a small system) and by checking of 
conserved quantities. 

Although we have here reported only on ID simulations, 
we have also implemented 2D and 3D calculations and 
found that the method works reliably in higher dimen- 
sions. Extensions of the method to implement s-wave scat- 
tering interactions will enable the study of non-equilibrium 
dynamics in a broader class of fermionic systems of cur- 
rent experimental interest, such as atomic Mott insulators 
in optical lattices and the BCS-BEC crossover problem. 

* * * 
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